********************************************************************************
*     
********************************************************************************
clear programs
*clear all
set more off

capture program drop main
program define main
	disp "Calling subprograms"
	*LoadData
	*EventStudyPlots
	RobustExcludeControls
	RobustHeatAltBaseline
	RobustTimeTrendsDiffer
	RobustClustering
	Reg_buildtype
	SaveCoefficients
	PreTrends
end

*===============================================================================
* 0) Load Data
capture program drop LoadData
program define LoadData
*===============================================================================
di "Load data"

cd "$GeneralModified"
use "BalancedData",replace

end

*===============================================================================
* 1) Main event study estimates
capture program drop EventStudyPlots
program define EventStudyPlots
*===============================================================================
di "EventStudyPlots"

eststo Zest2: quietly reg t_b3 t_b2_fm*,nocon
*eststo Zest3: quietly reg t_b4 t_b3_fm*,nocon

*===============================================================================
* Figure 1
*===============================================================================

eststo c1om_tps01_RR_1: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

*This uses balancedset_change, which exludes households after the reward value changes
eststo c1_tps01_RR_2: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges==1 & shortgapV2<=12 & balancedset_change==1 & tps==1 | (tps==0 & balancedset_change==1),fe vce(cl id)

eststo c2om_tps01_RR_3: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges>=2 & shortgapV2<=12 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

*c2 vs c3om, c3 vs c4om
*balancedset_change
eststo c2_tps01_RR_4: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges==2 & shortgapV2<=12 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

eststo c3om_tps01_RR_4: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges>=3 & shortgapV2<=12 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

eststo c3_tps01_RR_5: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges==3 & balancedset_change==1 & tps==1 & shortgapV2<=12 | (tps==0 & balancedset_change==1),fe vce(cl id)

eststo c4om_tps01_RR_5: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* reduc_gap_* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm* t_a6_fm* t_a7_fm* t_a8_fm* i.date if num_challenges>=4 &  balancedset_change==1 & tps==1 & shortgapV2<=12 | (tps==0 & balancedset_change==1),fe vce(cl id)

*esttab c1om_tps01_RR_1 c2om_tps01_RR_3,keep(t_b1_fm* t_a0_*) compress wide nogaps star(* 0.10 ** 0.05 *** 0.01) se order(t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12) 
*-------------------------------------------------------------------------------
*4 plots in one
*-------------------------------------------------------------------------------

*First plot
coefplot (c1om_tps01_RR_1,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) mcolor(navy) lc(gs4) msize(small) ciop(recast(rline) lc(gs4) lp(dash) msize(small))) ///
	(c1om_tps01_RR_1,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(line) mcolor(navy) lc(gs4) msize(small) ciop(recast(rline) lc(gs4) lp(dash) msize(small))) ///
	(c1om_tps01_RR_1,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small) alt) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xtick(12(12)120) xmtick(0(1)132) xlabel(none) ytitle("") xtitle("") legend(off) graphregion(color(white)) fysize(23.45) ///
	text(0.12 38 "(a) All households",size(medsmall)  justification(left) placement(e)) 
graph rename p1,replace	

*First plot—single version for appendix
coefplot (c1om_tps01_RR_1,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) mcolor(navy) lc(gs4) lw(medthick) ciop(recast(rline) lc(gs4) lp(dash) msize(small))) ///
	(c1om_tps01_RR_1,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(line) mcolor(navy) lc(gs4) lw(medthick) ciop(recast(rline) lc(gs4) lp(dash) msize(small))) ///
	(c1om_tps01_RR_1,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	xtick(12(12)120) legend(off) graphregion(color(white))
graph rename p1_appendix,replace	


*Second plot
coefplot (c1_tps01_RR_2,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c1_tps01_RR_2,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c1_tps01_RR_2,keep(t_a0_fm*) mc(maroon) lw(thick) recast(sc) msize(small) noci) ///
	(c2om_tps01_RR_3,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(forest_green) lp(dash) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash))) ///
	(c2om_tps01_RR_3,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) lp(dash) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash))) ///
	(c2om_tps01_RR_3,keep(t_a0_fm* t_a1_fm*) mc(maroon) msymbol(T) lw(thick) recast(sc) msize(small) noci) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("") xtitle("") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small) alt) ///
	xlabel(none) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	legend(off) fysize(23.45) text(0.125 38 "(b) Single vs two or more challenges",size(medsmall) justification(left) placement(e)) 
graph rename p2,replace

*Second plot, appendix
coefplot (c1_tps01_RR_2,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c1_tps01_RR_2,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c1_tps01_RR_2,keep(t_a0_fm*) mc(maroon) lw(thick) recast(sc) msize(small) noci) ///
	(c2om_tps01_RR_3,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lp(dash) lc(forest_green) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash))) ///
	(c2om_tps01_RR_3,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green)lp(dash) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash))) ///
	(c2om_tps01_RR_3,keep(t_a0_fm* t_a1_fm*) mc(maroon) lw(thick) recast(sc) msize(small) msymbol(T) noci) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin)) legend(order(2 "One Challenge" 7 "Two or more Challenges"))
graph rename p2_appendix,replace
cd "$SaveOutputFiles"
graph export p2_appendix.png,replace

*3rd plot
coefplot (c2_tps01_RR_4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) msize(small) lp(dot))) ///
	(c2_tps01_RR_4,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) lp(dot) msize(small))) ///
	(c2_tps01_RR_4,keep(t_a0_fm* t_a1_fm*) mc(maroon) lw(thick) recast(sc) noci msize(small)) ///
	(c3om_tps01_RR_4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c3om_tps01_RR_4,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c3om_tps01_RR_4,keep(t_a0_fm* t_a1_fm* t_a2_fm*) mc(maroon) msymbol(T) lw(thick) recast(sc) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("") xtitle("") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small) alt) ///
	xlabel(12(12)120 ) xmtick(0(1)132) xlabel(none) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	legend(off) title("") fysize(23.45) text(0.125 38 "(c) Two vs three or more challenges",size(medsmall) justification(left) placement(e)) 
graph rename p3,replace

*3rd plot, appendix
coefplot (c2_tps01_RR_4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) msize(small) lp(dot))) ///
	(c2_tps01_RR_4,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) lp(dot) msize(small))) ///
	(c2_tps01_RR_4,keep(t_a0_fm* t_a1_fm*) mc(maroon) lw(thick) recast(sc) noci msize(small)) ///
	(c3om_tps01_RR_4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lp(dash) lc(forest_green) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c3om_tps01_RR_4,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c3om_tps01_RR_4,keep(t_a0_fm* t_a1_fm* t_a2_fm*) mc(maroon) msymbol(T) lw(thick) recast(sc) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin)) legend(order(2 "Two Challenges" 7 "Three or more Challenges"))
graph rename p3_appendix,replace
cd "$SaveOutputFiles"
graph export p3_appendix.png,replace

*4th plot
coefplot (c3_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) lp(dot) msize(small))) ///
	(c3_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c3_tps01_RR_5,keep(t_a0_fm* t_a1_fm* t_a2_fm*) mc(maroon) lw(thick) recast(sc) noci) ///
	(c4om_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c4om_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c4om_tps01_RR_5,keep(t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm*) mc(maroon)msymbol(T) lw(thick) recast(sc) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small) alt) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin)) ///
	legend(off)  title("") fysize(29.65) text(0.125 38 "(d) Three vs four or more challenges",size(medsmall) justification(left) placement(e)) 
graph rename p4,replace

*4th plot-appendix
coefplot (c3_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) msize(small) ciop(recast(rline) lw(thick) lc(midblue) lp(dot) msize(small))) ///
	(c3_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) ciop(recast(rline) lw(thick) lc(midblue) lp(dot))) ///
	(c3_tps01_RR_5,keep(t_a0_fm* t_a1_fm* t_a2_fm*) mc(maroon) lw(thick) recast(sc) noci) ///
	(c4om_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c4om_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) lp(dash) msize(small) ciop(recast(rline) lw(medthick) lc(forest_green) lp(dash) msize(small))) ///
	(c4om_tps01_RR_5,keep(t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm*) mc(maroon) msymbol(T) lw(thick) recast(sc) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin)) legend(order(2 "Three Challenges" 7 "Four or more Challenges"))
graph rename p4_appendix,replace
cd "$SaveOutputFiles"
graph export p4_appendix.png,replace

graph combine p1 p2 p3 p4,r(4) graphregion(color(white)) iscale(0.6) imargin(vsmall) name(DID_RR_4plots_v2,replace)
graph display DID_RR_4plots_v2, xsize(8) ysize(9) margins(vsmall)
cd "$SaveOutputFiles"
graph export DID_RR_4plots_v2.png,replace

*------------------------------------------------------------------------------- 
*Plot all estimates overlapped
coefplot (c1_tps01_RR_2,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(midblue) noci) ///
	(c1_tps01_RR_2,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(midblue) noci) ///
	(c2_tps01_RR_4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(forest_green) noci) ///
	(c2_tps01_RR_4,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(forest_green) noci) ///
	(c3_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(red) noci) ///
	(c3_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(red) noci) ///
	(c4om_tps01_RR_5,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(magenta) noci) ///
	(c4om_tps01_RR_5,keep(t_b1_fm* t_a0_fm* t_a1_fm* t_a2_fm* t_a3_fm* t_a4_fm* t_a5_fm*) recast(line) lc(magenta) noci) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.20 0.1)) ylabel(-0.20 "-20%" -0.15 "-15%" -0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12)) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m*) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) ///
	legend(order(2 "One Challenge" 4 "Two Challenges" 6 "Three Challenges" 8 "Four Challenges"))
//OLD name: *COEFPLOT BALANCED SET R.13.20
graph rename SelectionAll_R_13_20_v2,replace
cd "$SaveOutputFiles"
graph export SelectionAll_R_13_20_v2.png,replace


end

*===============================================================================
* 2) Robustness—exclude non-participants & all non-participants
capture program drop RobustExcludeControls
program define RobustExcludeControls
*===============================================================================
di "RobustExcludeControls"
*Check robustness excluding control households
eststo c1om_tps1_RR_1: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM"),fe vce(cl id)

*First plot
coefplot (c1om_tps1_RR_1,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps1_RR_1,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps1_RR_1,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	xtick(12(12)120) legend(off) 
graph rename p1_tps1,replace
cd "$SaveOutputFiles"
graph export p1_tps1.png,replace

esttab c1om_tps01_RR_1 c1om_tps1_RR_1,keep(t_a0_*) compress wide nogaps star(* 0.10 ** 0.05 *** 0.01) se
//this shows that the estimates between these models are slightly different.

*-------------------------------------------------------------------------------
*Estimate using all non-participants

eststo c1om_tps01_RR_3: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0),fe vce(cl id)

*First plot
coefplot (c1om_tps01_RR_3,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_3,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_3,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	xtick(12(12)120) legend(off) graphregion(color(white))
graph rename p1_tps01_allcontrols,replace
cd "$SaveOutputFiles"
graph export p1_tps01_allcontrols.png,replace

end

*===============================================================================
* 3) Robustness—heating/altbaseline
capture program drop RobustHeatAltBaseline
program define RobustHeatAltBaseline
*===============================================================================
di "RobustHeatAltBaseline"
*Alternate baseline estimates
eststo c1om_tps01_RR_2: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b2_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & balancedset_all==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_all==1),fe vce(cl id)

eststo Zest3: quietly reg t_b4 t_b3_fm*,nocon

*Alternate baseline
coefplot (c1om_tps01_RR_2,keep(t_b5_fm* t_b4_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_2,keep(t_b2_fm* t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_2,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest3,keep(t_b3_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12) labsize(small)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	xtick(12(12)120) legend(off) graphregion(color(white))
graph rename p1_altbaseline,replace	
cd "$SaveOutputFiles"
graph export p1_altbaseline.png,replace


*-------------------------------------------------------------------------------
*Heating plots

*Updated heating 0 vs 1 plots

eststo c1om_tps01_RR_ht0: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==0 | (tps==0 & balancedset_change==1 & heating==0),fe vce(cl id)
eststo c1om_tps01_RR_ht1: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==1 | (tps==0 & balancedset_change==1 & heating==1),fe vce(cl id)

*Single heating plot before combining
coefplot (c1om_tps01_RR_ht0,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(navy) lw(medthick) ciop(recast(rline) lw(medthin) lc(navy) lp(dash))) ///
	(c1om_tps01_RR_ht0,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(line) lc(navy) lw(medthick) ciop(recast(rline) lw(medthin) lc(navy) lp(dash))) ///
	(c1om_tps01_RR_ht1,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(line) lc(maroon) lw(medthick) ciop(recast(rline) lw(thin) lc(maroon) lp(dash))) ///
	(c1om_tps01_RR_ht1,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(line) lc(maroon)  lw(medthick) ciop(recast(rline) lw(thin) lc(maroon) lp(dash))) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") graphregion(color(white)) ///
	yscale(range(-0.10 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12)) ///
	xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m*) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) ///
	legend(order(2 "Non-Electric Heating" 6 "Electric Space Heating") size(vsmall))
graph rename c1om_heat01,replace


*Plot the seasonal estimates
disp "Season effects"
gen monthofyear=month(dofm(date))

forvalues x=1/12 {
	gen t_a0_m`x'=t_a0 if monthofyear==`x'
	replace t_a0_m`x'=0 if t_a0_m`x'==.
}

*Use reference year of t_b2 for consistencey with the other heating estimates above 
forvalues i=1/12 {
	di `i'
	eststo c1om_sea_ht1_2_m`i': quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0_m`i' t_t_m* i.date if balancedset_change==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==1 & monthofyear==`i'| (tps==0 & balancedset_change==1 & heating==1 & monthofyear==`i'),fe vce(cl id)
}

 forvalues i=1/12 {
	di `i'
	eststo c1om_sea_ht0_2_m`i': quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0_m`i' t_t_m* i.date if balancedset_change==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==0 & monthofyear==`i'| (tps==0 & balancedset_change==1 & heating==0 & monthofyear==`i'),fe vce(cl id)
}
 
coefplot (c1om_sea_ht0_2_m1 \ c1om_sea_ht0_2_m2 \ c1om_sea_ht0_2_m3 \ c1om_sea_ht0_2_m4 \ c1om_sea_ht0_2_m5 \ c1om_sea_ht0_2_m6 \ c1om_sea_ht0_2_m7 \ c1om_sea_ht0_2_m8 \ c1om_sea_ht0_2_m9  \ c1om_sea_ht0_2_m10 \ c1om_sea_ht0_2_m11 \ c1om_sea_ht0_2_m12,keep(t_a0_m*) recast(con) ciop(recast(rline) lw(thick) lc(navy) lp(dot))) ///
	(c1om_sea_ht1_2_m1 \ c1om_sea_ht1_2_m2 \ c1om_sea_ht1_2_m3 \ c1om_sea_ht1_2_m4 \ c1om_sea_ht1_2_m5 \ c1om_sea_ht1_2_m6 \ c1om_sea_ht1_2_m7 \ c1om_sea_ht1_2_m8 \ c1om_sea_ht1_2_m9  \ c1om_sea_ht1_2_m10 \ c1om_sea_ht1_2_m11 \ c1om_sea_ht1_2_m12,keep(t_a0_m*) recast(con) ciop(recast(rline) lw(thick) lc(maroon) lp(dot))) ///
	,vert nooff yscale(range(-0.1 0)) ylabel(-0.1 "-10%" -0.08 "-8%" -0.06 "-6%" -0.04 "-4%" -0.02 "-2%" 0 "0%",glcolor(gs12)) ///
	ytitle("Change From Reference Year") graphregion(color(white)) ///
	legend(label(2 "Non-Electric Heating") label(4 "Electric Space Heating")) ///
	order(t_a0_m1 t_a0_m2 t_a0_m3 t_a0_m4 t_a0_m5 t_a0_m6 t_a0_m7 t_a0_m8 t_a0_m9 t_a0_m10 t_a0_m11 t_a0_m12) ///
	rename(t_a0_m1="JAN" t_a0_m2="FEB" t_a0_m3="MAR" t_a0_m4="APR" t_a0_m5="MAY" t_a0_m6="JUN" t_a0_m7="JUL" t_a0_m8="AUG" t_a0_m9="SEP" t_a0_m10="OCT" t_a0_m11="NOV" t_a0_m12="DEC")
*PLOT SEASONAL S.4.3
graph rename seasonal_EL_S4_3,replace

grc1leg seasonal_EL_S4_3 c1om_heat01,r(2) graphregion(color(white)) iscale(0.75) imargin(vsmall) name(SEASONAL_S_4_4,replace) legendfrom(c1om_heat01)
graph display SEASONAL_S_4_4, xsize(9) ysize(10)
cd "$SaveOutputFiles"
graph export SEASONAL_S_4_4.png,replace

*To compare heating breakdowns its helpful to compare the conservation vs. t_b1 not t_b2
forvalues i=1/12 {
	di `i'
	eststo c1om_sea_ht1_3_m`i': quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b2 t_a0_m`i' t_t_m* i.date if balancedset_change==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==1 & monthofyear==`i'| (tps==0 & balancedset_change==1 & heating==1 & monthofyear==`i'),fe vce(cl id)
}

 forvalues i=1/12 {
	di `i'
	eststo c1om_sea_ht0_3_m`i': quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b2 t_a0_m`i' t_t_m* i.date if balancedset_change==1 & begin_challenge_1<=monthly("2014m8","YM") & heating==0 & monthofyear==`i'| (tps==0 & balancedset_change==1 & heating==0 & monthofyear==`i'),fe vce(cl id)
}

coefplot (c1om_sea_ht0_3_m1 \ c1om_sea_ht0_3_m2 \ c1om_sea_ht0_3_m3 \ c1om_sea_ht0_3_m4 \ c1om_sea_ht0_3_m5 \ c1om_sea_ht0_3_m6 \ c1om_sea_ht0_3_m7 \ c1om_sea_ht0_3_m8 \ c1om_sea_ht0_3_m9  \ c1om_sea_ht0_3_m10 \ c1om_sea_ht0_3_m11 \ c1om_sea_ht0_3_m12,keep(t_a0_m*) recast(con) ciop(recast(rline) lw(thick) lc(navy) lp(dot))) ///
	(c1om_sea_ht1_3_m1 \ c1om_sea_ht1_3_m2 \ c1om_sea_ht1_3_m3 \ c1om_sea_ht1_3_m4 \ c1om_sea_ht1_3_m5 \ c1om_sea_ht1_3_m6 \ c1om_sea_ht1_3_m7 \ c1om_sea_ht1_3_m8 \ c1om_sea_ht1_3_m9  \ c1om_sea_ht1_3_m10 \ c1om_sea_ht1_3_m11 \ c1om_sea_ht1_3_m12,keep(t_a0_m*) recast(con) ciop(recast(rline) lw(thick) lc(maroon) lp(dot))) ///
	,vert nooff yscale(range(-0.1 0)) ylabel(-0.1 "-10%" -0.08 "-8%" -0.06 "-6%" -0.04 "-4%" -0.02 "-2%" 0 "0%",glcolor(gs12)) ///
	ytitle("Change From Reference Year") graphregion(color(white)) ///
	legend(label(2 "Non-Electric Heating") label(4 "Electric Space Heating")) ///
	order(t_a0_m1 t_a0_m2 t_a0_m3 t_a0_m4 t_a0_m5 t_a0_m6 t_a0_m7 t_a0_m8 t_a0_m9 t_a0_m10 t_a0_m11 t_a0_m12) ///
	rename(t_a0_m1="JAN" t_a0_m2="FEB" t_a0_m3="MAR" t_a0_m4="APR" t_a0_m5="MAY" t_a0_m6="JUN" t_a0_m7="JUL" t_a0_m8="AUG" t_a0_m9="SEP" t_a0_m10="OCT" t_a0_m11="NOV" t_a0_m12="DEC")
graph rename seasonal_EL_S4_3,replace

*esttab c1om_sea_ht0*,keep(t_a0_m*)compress nogaps star(* 0.10 ** 0.05 *** 0.01) se
*esttab c1om_sea_ht1*,keep(t_a0_m*)compress nogaps star(* 0.10 ** 0.05 *** 0.01) se
*esttab c1om_sea_ht1_2_m6 c1om_sea_ht0_2_m6,keep(t_*)compress wide nogaps star(* 0.10 ** 0.05 *** 0.01) se
*esttab c1om_sea_ht1_2_m1 c1om_sea_ht0_2_m1,keep(t_*)compress wide nogaps star(* 0.10 ** 0.05 *** 0.01) se

*Calculate the household average use (pre-program) over the 12 months of the year
gen month=month(dofm(date))
gen year=year(dofm(date))
tabstat kwhw if year==2006 & tps==1 & heating==0 & balancedset_change==1,by(month)
tabstat kwhw if year==2006 & tps==1 & heating==1 & balancedset_change==1,by(month)

end

*===============================================================================
* 4) Robustness—Time trends
capture program drop RobustTimeTrendsDiffer
program define RobustTimeTrendsDiffer
*===============================================================================
*Define date as separate trends depending on characteristics like building type and heating
di "RobustTimeTrendsDiffer"
cd "$GeneralModified"
use "BalancedData",replace

eststo Zest2: quietly reg t_b3 t_b2_fm*,nocon

br id tps date buildingtype heating
egen date_ht=group(date heating)

eststo c1om_tps01_RR_TT_1: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date_ht if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

coefplot (c1om_tps01_RR_TT_1,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_TT_1,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_TT_1,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) ///
		xlabel(12(12)120 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	ytitle("Change From Reference Year") xtitle("") legend(off) graphregion(color(white)) title("Trend by date-heating",size(medsmall)) 
graph rename c1om_tps01_RR_TT_1,replace
cd "$SaveOutputFiles"
graph export c1om_tps01_RR_TT_1.png,replace

*heating and building type separate trends
egen date_htbt=group(date heating buildingtype)

eststo c1om_tps01_RR_TT_2: quietly xtreg lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* i.date_htbt if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),fe vce(cl id)

coefplot (c1om_tps01_RR_TT_2,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_TT_2,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_TT_2,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) ///
	xlabel(12(12)120 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	ytitle("Change From Reference Year") xtitle("") legend(off) graphregion(color(white)) title("Trend by date-heating-building type",size(medsmall)) 
graph rename c1om_tps01_RR_TT_2,replace
cd "$SaveOutputFiles"
graph export c1om_tps01_RR_TT_2.png,replace

end

*===============================================================================
* 5) Robustness—twoway clustering
capture program drop RobustClustering
program define RobustClustering
*===============================================================================
di "RobustClustering"
cd "$GeneralModified"
use "BalancedData",replace

eststo Zest2: quietly reg t_b3 t_b2_fm*,nocon

eststo c1om_tps01_RR_C4: quietly reghdfe lnkwhw t_b10_fm* t_b9_fm* t_b8_fm* t_b7_fm* t_b6_fm* t_b5_fm* t_b4_fm* t_b3_fm* t_b1_fm* t_a0_fm* t_t_m* if num_challenges>=1 & balancedset_change==1 & tps==1 & begin_challenge_1<=monthly("2014m8","YM") | (tps==0 & balancedset_change==1),absorb(id date) vce(cl i.id i.date)

coefplot (c1om_tps01_RR_C4,keep(t_b5_fm* t_b4_fm* t_b3_fm*) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_C4,keep(t_b1_fm* t_a0_fm* t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) recast(con) mcolor(navy) lc(gs4) msize(small) ciop(recast(rcap) lc(gs4) lp(longdash) msize(small))) ///
	(c1om_tps01_RR_C4,keep(t_a0_fm*) recast(sc) mc(maroon) noci msize(small)) ///
	(Zest2,keep(t_b2_fm*) recast(sc) msymbol(i) noci) ///
	,vert nooff ///
	yscale(range(-0.1 0.1)) ylabel(-0.1 "-10%" -0.05 "-5%" 0 "0%" 0.05 "5%" 0.1 "10%",glcolor(gs12)) ///
	order(t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) ///
	xline(12,lp(dot)) xline(24,lp(dot)) xline(36,lp(dot)) xline(48,lp(dot)) xline(60,lp(dot)) xline(72,lp(dot)) xline(84,lp(dot)) xline(96,lp(dot)) xline(108,lp(dot)) xline(120,lp(dot)) yline(0,lcolor(black) lwidth(medthin))  ///
	xtick(12(12)120) xmtick(0(1)132) xlabel(12(12)36 12 `"-48"' 24 `"-36"' 36 `"-24"' 48 `"-12"' 60 `"Start"' 72 `"12"' 84 `"24"' 96 `"36"' 108 `"48"' 120 `"60"' ) xmtick(0(1)132) ///
	ytitle("Change From Reference Year") xtitle("Months Pre/Post 1st Challenge Start") legend(off) graphregion(color(white)) title("")
graph rename p1_doublecluster4,replace
cd "$SaveOutputFiles"
graph export p1_doublecluster4.png,replace

end

*===============================================================================
* Event study annual specification for household characteristics
capture program drop Reg_buildtype
program define Reg_buildtype
*===============================================================================
di "Reg_buildtype"
cd "$GeneralModified"
use "BalancedData",replace

*Simplify dataset; estimate annual specification but using monthly data
keep id tps date uniid reduc kwhw lnkwhw num_challenges t_a0 t_b* t_t_m* buildingtype heating begin_challenge_1 balancedset_change

*Merge in info needed like total_kwh_num
cd "$datafrom"
merge m:1 id using unique_buildingchars_class,keepusing(floorarea value_tot) keep(1 3) nogen

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Estimates by building type
tab buildingtype,nolabel
local char="buildingtype"
eststo annual_tps01_char_m1_1: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==1,fe vce(cl id)
eststo annual_tps01_char_m1_2: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==2,fe vce(cl id)
eststo annual_tps01_char_m1_3: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==3,fe vce(cl id)
eststo annual_tps01_char_m1_4: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==4,fe vce(cl id)
eststo annual_tps01_char_m1_5: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==5,fe vce(cl id)
eststo annual_tps01_char_m1_999: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==999,fe vce(cl id)

estadd local buildtype "1-SFD" , replace : annual_tps01_char_m1_1
estadd local buildtype "2-SFD" , replace : annual_tps01_char_m1_2
estadd local buildtype "Apartment" , replace : annual_tps01_char_m1_4
estadd local buildtype "Townhouse" , replace : annual_tps01_char_m1_5

esttab annual_tps01_char_m1_*,keep(t_a0) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se
esttab annual_tps01_char_m1_*,keep(t_a0) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se noobs tex fragment s(buildtype)

*Number of participants and non-participants by hosuehold type
tabstat id if uniid==1 & balancedset_change==1 & tps==1,by(buildingtype) s(N)
tabstat id if uniid==1 & balancedset_change==1 & tps==0,by(buildingtype) s(N)

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*Estimates by quartiles of pre-program use

*Generate xtiles based on 2006
gen year=year(dofm(date))
capture drop temp
bysort id: egen temp=mean(kwhw) if year==2006
bysort id: egen kwhw2006=mean(temp)

capture drop temp
xtile quant_kwhw2006=kwhw2006 if balancedset_change==1 & tps==0,n(4)
xtile temp=kwhw2006 if balancedset_change==1 & tps==1,n(4)
replace quant_kwhw2006=temp if tps==1

local char="quant_kwhw2006"
eststo annual_tps01_char_m4_1: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==1,fe vce(cl id)
eststo annual_tps01_char_m4_2: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==2,fe vce(cl id)
eststo annual_tps01_char_m4_3: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==3,fe vce(cl id)
eststo annual_tps01_char_m4_4: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==4,fe vce(cl id)

esttab annual_tps01_char_m4_*,keep(t_a0) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se

tabstat id if uniid==1 & balancedset_change==1 & tps==1,by(quant_kwhw2006) s(N)
tabstat id if uniid==1 & balancedset_change==1 & tps==0,by(quant_kwhw2006) s(N)

tabstat kwhw if uniid==1 & balancedset_change==1 & tps==1 & year==2006,by(quant_kwhw2006)

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*quartiles of assessed value
capture drop temp
xtile quart_value=value_tot if balancedset_change==1 & tps==0,n(4)
xtile temp=value_tot if balancedset_change==1 & tps==1,n(4)
replace quart_value=temp if tps==1

local char="quart_value"
eststo annual_tps01_char_m2_1: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==1,fe vce(cl id)
eststo annual_tps01_char_m2_2: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==2,fe vce(cl id)
eststo annual_tps01_char_m2_3: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==3,fe vce(cl id)
eststo annual_tps01_char_m2_4: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==4,fe vce(cl id)

esttab annual_tps01_char_m2_*,keep(t_a0) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se
tabstat tps if uniid==1 & balancedset_change==1 & tps==1,by(quart_value) s(N)
tabstat value_tot if uniid==1 & balancedset_change==1 & tps==1,by(quart_value)

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
*quartlies of floor area 
capture drop temp
xtile quart_floorarea=floorarea if balancedset_change==1 & tps==0,n(4)
xtile temp=floorarea if balancedset_change==1 & tps==1,n(4)
replace quart_floorarea=temp if tps==1

local char="quart_floorarea"
eststo annual_tps01_char_m3_1: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==1,fe vce(cl id)
eststo annual_tps01_char_m3_2: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==2,fe vce(cl id)
eststo annual_tps01_char_m3_3: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==3,fe vce(cl id)
eststo annual_tps01_char_m3_4: quietly xtreg lnkwhw t_b10 t_b9 t_b8 t_b7 t_b6 t_b5 t_b4 t_b3 t_b1 t_a0 t_t_m* i.date if balancedset_change==1 & `char'==4,fe vce(cl id)

esttab annual_tps01_char_m3_*,keep(t_a0) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se
tabstat floorarea if uniid==1 & balancedset_change==1 & tps==1,by(quart_floorarea)

end

*===============================================================================
* Output of table of coefficient estimates for the appendix
capture program drop SaveCoefficients
program define SaveCoefficients

di "SaveCoefficients"

*Save estimates
/*
c1om_tps01_RR_1
c1om_tps01_RR_3
c1om_tps01_RR_ht0
c1om_tps01_RR_ht1

c1_tps01_RR_2
c2om_tps01_RR_3
c2_tps01_RR_4
c3om_tps01_RR_4
c3_tps01_RR_5
c4om_tps01_RR_5
*/

cd "$SaveOutputFiles"

esttab c1om_tps01_RR_1 c1om_tps01_RR_3 c1om_tps01_RR_ht0 c1om_tps01_RR_ht1 using c1om_heat01.tex,drop(*.date) compress se nogaps mti order(t_b10_fm1 t_b10_fm2 t_b10_fm3 t_b10_fm4 t_b10_fm5 t_b10_fm6 t_b10_fm7 t_b10_fm8 t_b10_fm9 t_b10_fm10 t_b10_fm11 t_b10_fm12 t_b9_fm1 t_b9_fm2 t_b9_fm3 t_b9_fm4 t_b9_fm5 t_b9_fm6 t_b9_fm7 t_b9_fm8 t_b9_fm9 t_b9_fm10 t_b9_fm11 t_b9_fm12 t_b8_fm1 t_b8_fm2 t_b8_fm3 t_b8_fm4 t_b8_fm5 t_b8_fm6 t_b8_fm7 t_b8_fm8 t_b8_fm9 t_b8_fm10 t_b8_fm11 t_b8_fm12 t_b7_fm1 t_b7_fm2 t_b7_fm3 t_b7_fm4 t_b7_fm5 t_b7_fm6 t_b7_fm7 t_b7_fm8 t_b7_fm9 t_b7_fm10 t_b7_fm11 t_b7_fm12 t_b6_fm1 t_b6_fm2 t_b6_fm3 t_b6_fm4 t_b6_fm5 t_b6_fm6 t_b6_fm7 t_b6_fm8 t_b6_fm9 t_b6_fm10 t_b6_fm11 t_b6_fm12 t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 t_t_m1 t_t_m2 t_t_m3 t_t_m4 t_t_m5 t_t_m6 t_t_m7 t_t_m8 t_t_m9 t_t_m10 t_t_m11 t_t_m12 t_t_m13 t_t_m14 t_t_m15 t_t_m16 t_t_m17 t_t_m18 t_t_m19 t_t_m20 t_t_m21 t_t_m22 t_t_m23 t_t_m24 t_t_m25 t_t_m26 t_t_m27 t_t_m28 t_t_m29 t_t_m30 t_t_m31 t_t_m32 t_t_m33 t_t_m34 t_t_m35 t_t_m36 t_t_m37 t_t_m38 t_t_m39 t_t_m40 t_t_m41 t_t_m42 t_t_m43 t_t_m44 t_t_m45 t_t_m46 t_t_m47 t_t_m48 t_t_m49 t_t_m50 t_t_m51 t_t_m52 t_t_m53 t_t_m54 t_t_m55 t_t_m56 t_t_m57 t_t_m58 t_t_m59 t_t_m60) tex replace


esttab c1_tps01_RR_2 c2om_tps01_RR_3 c2_tps01_RR_4 c3om_tps01_RR_4 c3_tps01_RR_5 c4om_tps01_RR_5 using c1om.tex,drop(*.date) b(a2) compress se nogaps mti order(t_b10_fm1 t_b10_fm2 t_b10_fm3 t_b10_fm4 t_b10_fm5 t_b10_fm6 t_b10_fm7 t_b10_fm8 t_b10_fm9 t_b10_fm10 t_b10_fm11 t_b10_fm12 t_b9_fm1 t_b9_fm2 t_b9_fm3 t_b9_fm4 t_b9_fm5 t_b9_fm6 t_b9_fm7 t_b9_fm8 t_b9_fm9 t_b9_fm10 t_b9_fm11 t_b9_fm12 t_b8_fm1 t_b8_fm2 t_b8_fm3 t_b8_fm4 t_b8_fm5 t_b8_fm6 t_b8_fm7 t_b8_fm8 t_b8_fm9 t_b8_fm10 t_b8_fm11 t_b8_fm12 t_b7_fm1 t_b7_fm2 t_b7_fm3 t_b7_fm4 t_b7_fm5 t_b7_fm6 t_b7_fm7 t_b7_fm8 t_b7_fm9 t_b7_fm10 t_b7_fm11 t_b7_fm12 t_b6_fm1 t_b6_fm2 t_b6_fm3 t_b6_fm4 t_b6_fm5 t_b6_fm6 t_b6_fm7 t_b6_fm8 t_b6_fm9 t_b6_fm10 t_b6_fm11 t_b6_fm12 t_b5_fm1 t_b5_fm2 t_b5_fm3 t_b5_fm4 t_b5_fm5 t_b5_fm6 t_b5_fm7 t_b5_fm8 t_b5_fm9 t_b5_fm10 t_b5_fm11 t_b5_fm12 t_b4_fm1 t_b4_fm2 t_b4_fm3 t_b4_fm4 t_b4_fm5 t_b4_fm6 t_b4_fm7 t_b4_fm8 t_b4_fm9 t_b4_fm10 t_b4_fm11 t_b4_fm12 t_b3_fm1 t_b3_fm2 t_b3_fm3 t_b3_fm4 t_b3_fm5 t_b3_fm6 t_b3_fm7 t_b3_fm8 t_b3_fm9 t_b3_fm10 t_b3_fm11 t_b3_fm12 t_b2_fm1 t_b2_fm2 t_b2_fm3 t_b2_fm4 t_b2_fm5 t_b2_fm6 t_b2_fm7 t_b2_fm8 t_b2_fm9 t_b2_fm10 t_b2_fm11 t_b2_fm12 t_b1_fm1 t_b1_fm2 t_b1_fm3 t_b1_fm4 t_b1_fm5 t_b1_fm6 t_b1_fm7 t_b1_fm8 t_b1_fm9 t_b1_fm10 t_b1_fm11 t_b1_fm12 t_a0_fm1 t_a0_fm2 t_a0_fm3 t_a0_fm4 t_a0_fm5 t_a0_fm6 t_a0_fm7 t_a0_fm8 t_a0_fm9 t_a0_fm10 t_a0_fm11 t_a0_fm12 reduc_gap_1 t_a1_fm1 t_a1_fm2 t_a1_fm3 t_a1_fm4 t_a1_fm5 t_a1_fm6 t_a1_fm7 t_a1_fm8 t_a1_fm9 t_a1_fm10 t_a1_fm11 t_a1_fm12 reduc_gap_2 t_a2_fm1 t_a2_fm2 t_a2_fm3 t_a2_fm4 t_a2_fm5 t_a2_fm6 t_a2_fm7 t_a2_fm8 t_a2_fm9 t_a2_fm10 t_a2_fm11 t_a2_fm12 reduc_gap_3 t_a3_fm1 t_a3_fm2 t_a3_fm3 t_a3_fm4 t_a3_fm5 t_a3_fm6 t_a3_fm7 t_a3_fm8 t_a3_fm9 t_a3_fm10 t_a3_fm11 t_a3_fm12 reduc_gap_4 t_a4_fm1 t_a4_fm2 t_a4_fm3 t_a4_fm4 t_a4_fm5 t_a4_fm6 t_a4_fm7 t_a4_fm8 t_a4_fm9 t_a4_fm10 t_a4_fm11 t_a4_fm12 reduc_gap_5 t_a5_fm1 t_a5_fm2 t_a5_fm3 t_a5_fm4 t_a5_fm5 t_a5_fm6 t_a5_fm7 t_a5_fm8 t_a5_fm9 t_a5_fm10 t_a5_fm11 t_a5_fm12 reduc_gap_6 t_a6_fm1 t_a6_fm2 t_a6_fm3 t_a6_fm4 t_a6_fm5 t_a6_fm6 t_a6_fm7 t_a6_fm8 t_a6_fm9 t_a6_fm10 t_a6_fm11 t_a6_fm12 reduc_gap_7 t_a7_fm1 t_a7_fm2 t_a7_fm3 t_a7_fm4 t_a7_fm5 t_a7_fm6 t_a7_fm7 t_a7_fm8 t_a7_fm9 t_a7_fm10 t_a7_fm11 t_a7_fm12 reduc_gap_8 t_a8_fm1 t_a8_fm2 t_a8_fm3 t_a8_fm4 t_a8_fm5 t_a8_fm6 t_a8_fm7 t_a8_fm8 t_a8_fm9 t_a8_fm10 t_a8_fm11 t_a8_fm12) tex replace

end

*===============================================================================
* 1) Pre-treatment trends
capture program drop PreTrends
program define PreTrends
*===============================================================================
di "Pretrends"
cd "$GeneralModified"
use "BalancedData",replace

gen date_center=date-monthly("2006m1","YM")
gen TPSdate_center=date_center if tps==1
replace TPSdate_center=0 if TPSdate_center==.

*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*Pre trends, but using balancedset_change
eststo trend_ln_tps01_m8: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_change==1 & (begin_challenge_1>monthly("2008m6","YM") | tps==0) & date<monthly("2008m1","YM"),nocon vce(cl id)
eststo trend_ln_tps01_m9: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_change==1 & (begin_challenge_1>monthly("2009m6","YM") | tps==0) & date<monthly("2009m1","YM"),nocon vce(cl id)
eststo trend_ln_tps01_m10: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_change==1 & (begin_challenge_1>monthly("2010m6","YM") | tps==0) & date<monthly("2010m1","YM"),nocon vce(cl id)
eststo trend_ln_tps01_m11: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_change==1 & (begin_challenge_1>monthly("2011m6","YM") | tps==0) & date<monthly("2011m1","YM"),nocon vce(cl id)
eststo trend_ln_tps01_m12: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_change==1 & (begin_challenge_1>monthly("2012m6","YM") | tps==0) & date<monthly("2012m1","YM"),nocon vce(cl id)

*eststo trend_ln_tps01_m8: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_all==1 & (begin_challenge_1<=monthly("2014m8","YM") | tps==0) & (begin_challenge_1>monthly("2008m6","YM") | tps==0) & date<monthly("2008m1","YM"),nocon vce(cl id)
*eststo trend_ln_tps01_m9: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_all==1 & (begin_challenge_1<=monthly("2014m8","YM") | tps==0) & (begin_challenge_1>monthly("2009m6","YM") | tps==0) & date<monthly("2009m1","YM"),nocon vce(cl id)
*eststo trend_ln_tps01_m10: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_all==1 & (begin_challenge_1<=monthly("2014m8","YM") | tps==0) & (begin_challenge_1>monthly("2010m6","YM") | tps==0) & date<monthly("2010m1","YM"),nocon vce(cl id)
*eststo trend_ln_tps01_m11: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_all==1 & (begin_challenge_1<=monthly("2014m8","YM") | tps==0) & (begin_challenge_1>monthly("2011m6","YM") | tps==0) & date<monthly("2011m1","YM"),nocon vce(cl id)
*eststo trend_ln_tps01_m12: quietly reg lnkwhw date_center TPSdate_center ibn.tps if balancedset_all==1 & (begin_challenge_1<=monthly("2014m8","YM") | tps==0) & (begin_challenge_1>monthly("2012m6","YM") | tps==0) & date<monthly("2012m1","YM"),nocon vce(cl id)

estadd local before "2008" , replace: trend_ln_tps01_m8
estadd local before "2009" , replace: trend_ln_tps01_m9
estadd local before "2010" , replace: trend_ln_tps01_m10
estadd local before "2011" , replace: trend_ln_tps01_m11
estadd local before "2012" , replace: trend_ln_tps01_m12
estadd local balanced "Change" , replace: trend_ln_tps01_m*

*Add the total number of accounts
count if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2008m6","YM") | tps==0)
estadd scalar numHH=r(N), replace: trend_ln_tps01_m8

count if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2009m6","YM") | tps==0)
estadd scalar numHH=r(N), replace: trend_ln_tps01_m9

count if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2010m6","YM") | tps==0)
estadd scalar numHH=r(N), replace: trend_ln_tps01_m10

count if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2011m6","YM") | tps==0)
estadd scalar numHH=r(N), replace: trend_ln_tps01_m11

count if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2012m6","YM") | tps==0)
estadd scalar numHH=r(N), replace: trend_ln_tps01_m12

tab tps if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2008m6","YM") | tps==0),matcell(testmat)
estadd scalar numtps0=testmat[1,1], replace: trend_ln_tps01_m8
estadd scalar numtps1=testmat[2,1], replace: trend_ln_tps01_m8

tab tps if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2009m6","YM") | tps==0),matcell(testmat)
estadd scalar numtps0=testmat[1,1], replace: trend_ln_tps01_m9
estadd scalar numtps1=testmat[2,1], replace: trend_ln_tps01_m9

tab tps if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2010m6","YM") | tps==0),matcell(testmat)
estadd scalar numtps0=testmat[1,1], replace: trend_ln_tps01_m10
estadd scalar numtps1=testmat[2,1], replace: trend_ln_tps01_m10

tab tps if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2011m6","YM") | tps==0),matcell(testmat)
estadd scalar numtps0=testmat[1,1], replace: trend_ln_tps01_m11
estadd scalar numtps1=testmat[2,1], replace: trend_ln_tps01_m11

tab tps if uniid==1 & balancedset_change==1 & (begin_challenge_1>monthly("2012m6","YM") | tps==0),matcell(testmat)
estadd scalar numtps0=testmat[1,1], replace: trend_ln_tps01_m12
estadd scalar numtps1=testmat[2,1], replace: trend_ln_tps01_m12

estadd local years "2" , replace: trend_ln_tps01_m8
estadd local years "3" , replace: trend_ln_tps01_m9
estadd local years "4" , replace: trend_ln_tps01_m10
estadd local years "5" , replace: trend_ln_tps01_m11
estadd local years "6" , replace: trend_ln_tps01_m12
*REGRESSION_DID_SUMMARYSTATS T.6.3
esttab trend_ln_tps01_m8 trend_ln_tps01_m9 trend_ln_tps01_m10 trend_ln_tps01_m11 trend_ln_tps01_m12, compress nogaps star(* 0.10 ** 0.05 *** 0.01) se s(before numHH numtps0 numtps1 N, label("Before" "Households" "tps0" "tps1"))

esttab trend_ln_tps01_m8 trend_ln_tps01_m9 trend_ln_tps01_m10 trend_ln_tps01_m11 trend_ln_tps01_m12, b(%8.5f) drop(*.tps) compress nogaps star(* 0.10 ** 0.05 *** 0.01) se s(years numtps0 numtps1 N, fmt(%12.0f) label("Years" "Non-Participant IDs" "Participant IDs" "Observations")) tex fragment

end


*========================================================================
*Call the main program to run this file.
*========================================================================
main
